Microscopic origins of the anomalous melting behaviour of high-pressure sodium 
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Recent experiments have shown that sodium, 
a prototype simple metal at ambient conditions, 
exhibits unexpected complexity under high pres- 
surei— . One of the most puzzling phenomena in 
the behaviour of dense sodium is the pressure- 
induced drop in its melting temperature, which 
extends from 1000 K at ^30 GPa to as low as 
room temperature at ~^120 GPa^. Despite signif- 
icant theoretical effort to understand the anoma- 
lous melting'^-^- its origins have remained unclear. 
In this work, we reconstruct the sodium phase di- 
agram using an a6-mitio-quality neural-network 
potential-^'lS. We demonstrate that the reentrant 
behaviour results from the screening of interi- 
onic interactions by conduction electrons, which 
at high pressure induces a softening in the short- 
range repulsion. It is expected that such an ef- 
fect plays an important role in governing the be- 
haviour of a wide range of metals and alloys. 

At ambient conditions, sodium crystallizes in the body- 
centered cubic (bcc) structure, which transforms to the 
face-centered cubic (fee) phase upon compression to 
65 GPa^. At 105 GPa, the fee phase undergoes a trans- 
formation to the cI16 phase - a distorted variant of the 
bcc structure^. Further compression results in the forma- 
tion of a large variety of complex phases^, the existence 
of which is quite intriguing. However, the focus of this 
work is one of the most striking features of the sodium 
phase diagram - the anomalous melting curve. Diffrac- 
tion measurements have revealed a maximum in melt- 
ing temperature around 1000 K at ^30 GPa followed by 
a pressure-induced drop, which extends to nearly room 
temperature at '^120 GPai and spans the regions of sta- 
bility of three solid phases (Fig. [1^). This anomaly im- 
plies that in this range the liquid is denser than the solid. 

Although previous ah initio studies^^— have repro- 
duced the reentrant melting behaviour of sodium the 
physical origin of this anomaly still remains unclear. The 
hypothesis of structural and electronic transitions in liq- 
uid sodium proposed in one of the studies^ has not been 
confirmed by other authors^. The discrepancy between 
these simulation results has been attributed to a different 
treatment of the semicore states^. 

It is important to point out that due to the high com- 
putational cost of ah initio MD all simulations so far re- 
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FIG. 1. Sodium phase diagram, a. The squares represent 
experimentally measured melting points. The lines show ten- 
tative phase boundaries^, b. Theoretical results based on the 
NN potential. The squares represent coexistence points com- 
puted with thermodynamic integration. The lines show coex- 
istence curves traced by integrating the Clausius-Clapeyron 
equation (see Methods). 



ported^"— have used small simulation cells and estimated 
melting temperatures (Tm) by heating the solid phase un- 
til it melts. In such simulations, the absence of nucleation 
centers and fast heating schedules delay melting relative 
to the thermodynamic coexistence point and, therefore, 
only an upper bound on can be obtained. A proper re- 
construction of the melting curve requires comparing the 
free energies of the liquid and solid phases and involves 
long (^10 ns) MD simulations on large systems (~ 10'^ 
atoms). Performing such simulations with direct ah ini- 
tio methods is computationally too expensive at present. 
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Although several authors^iiS have shown that simplified 
models can also reproduce the anomalous melting be- 
havior the relation of these models to the actual physical 
interactions in sodium has remained unclear. 

In a previous worki^. we have shown that a neural- 
network (NN) potential for sodium based on well- 
converged electronic structure calculations, retains the 
accuracy of ab initio simulations at a highly reduced com- 
putational cost. The ability of the NN potential to re- 
produce numerous ab initio and experimental properties 
of the solid and the liquid phases in the pressure range 
up to 140 GPa has been shown in Ref. [l^- Addition- 
ally, Supplementary Fig. 1 demonstrates high accuracy 
of the NN potential in reproducing the ab initio energies 
for liquid sodium structures along MD trajectories in the 
relevant P-T range. 

Here wc utilize this NN potential to construct the 
sodium phase diagram. The phase diagram shown in 
Fig. [T|d is calculated using thermodynamic integration 
and by tracing the coexistence curve using the Clausius- 
Clapeyron equation (see Methods). The NN potential 
captures the experimentally observed sequence of the 
solid-state phase transitions bcc^fcc^-cI16 and the re- 
gions of stability of each phase. The shape of the melting 
curve with the maximum around 30 GPa is also correctly 
reproduced. The comparison of NN and DFT transition 
pressures at zero temperature (Fig. 6 in Ref. [l^ as well 
as the "heat- until- it-melts " melting curves obtained from 
NN and ab initio simulations (Supplementary Fig. 2) 
show that the quantitative errors in coexistence lines can 
be attributed to inaccuracies of the density functional 
and not to the errors of the NN fitting. Supplemen- 
tary Fig. 2 also demonstrates that Tm obtained with a 
non-equilibrium approach is overestimated compared to 
the proper thermodynamic calculation. To the best of 
our knowledge this is the first phase diagram calculated 
from a high-quality ab initio method. A careful analysis 
of structural and electronic properties of liquid sodium 
shows no evidence of liquid-liquid phase transitions pro- 
posed earlier— (see Supplementary Information). 

While the NN simulations reproduce adequately the 
sodium melting curve the physical origins of its reentrant 
behaviour are difficult to extract. To this effect, wc con- 
structed a density dependent pair potential based on the 
jellium model for the conducting electrons. According to 
this model the ionic structure of a metal can be replaced 
to a first approximation, by a uniform positively-charged 
background while conduction electrons can be treated as 
a uniform electron gas. The granularity of sodium is then 
introduced at the level of two-body interactions by im- 
mersing two Na atoms in such a jellium and calculating 
the energy as a function of the interatomic distance (see 
Methods). The effect of the pressure is taken into account 
by varying the density of the electron gas pe according 
to the NN equation of state. By doing so, we retrace a 
time-honoured idea used to construct effective potentials 
describing screening effects in metals by linear response 
theoryi^. Here we use a non-linear approach since it is 
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FIG. 2. The "heat-until- it-melts" curves obtained with the 
NN potential (green), effective pair potential based on the 
jellium model (red) and repulsive wall of the effective pair 
potential (blue). The points below 90 GPa are obtained by 
melting the bcc phase, the points above 90 GPa by melting 
the fee phase. 

expected that non-linear screening effects become impor- 
tant as the ions are brought closer to each other by the 
applied pressure. 

We found that the two-body effective potential ob- 
tained from the jellium model is well reproduced by a 
sum of the repulsive Yukawa potential and the oscilla- 
tory Friedcl term (Supplementary Fig. 3)^^^: 

Aexp(-fcor) Bcos[2fcj-(r - To)] 
'^W = + ^ ' (1) 

where kp is the Fermi wavenumber and A, B, ko, tq and 
m are density dependent fitting parameters (see Meth- 
ods). This analytical form of the potential is inspired 
by a linear screening theorjJ^ and reflects the nature of 
physical interactions in metallic sodium, in which the di- 
rect Na-Na repulsion is screened by the free-electron gas. 
In particular, the oscillatory term has its origin in the 
sharpness of the Fermi surface. 

MD simulations based on the effective pair potential re- 
produce semi-quantitatively the maximum in the melting 
curve (Fig. [2]) • The potential is also remarkably accurate 
in predicting the radial distribution functions (RDFs) of 
the liquid for pressures up to 100 GPa (Supplementary 
Fig. 4) . We thus believe that this effective density depen- 
dent two-body potential is able to capture the physics of 
the problem and can be used to shed light on the origin 
of the anomalous melting behaviour. It is worth men- 
tioning that the validity of this potential is limited to 
the region of pressures below about 120 GPa because be- 
yond this region the electron are promoted to the d-states 
and the simple free-electron two-body model is no longer 
appropriate. 
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We have examined the effect of two different parts of 
the effective potential on the melting curve. First, we 
considered the Yukawa part of the potential and, second, 
the repulsive wall of The latter was constructed by 
truncating the total potential at its first minimum and 
shifting it so that both the energy and force are zero at 
the truncation pointiii. We found that the Yukawa poten- 
tial cannot reproduce the reentrant melting behaviour. 
On the contrary, the melting curve obtained with the re- 
pulsive wall exhibits a maximum similar to that obtained 
with the full potential (Fig. [2]). This observation suggests 
that it is the short range repulsive part of the potential 
that provides the origin of the anomalous melting (see 
Supplementary information for a detailed pressure de- 
composition analysis that confirms this conclusion). 

The log-log plot of the interatomic force in the region of 
the repulsive potential wall shown in Fig. [3]) supports this 
point of view. It is seen that there is a large change of the 
slope as the potential changes from 1 / r^^ to a much softer 

behaviour (i.e. the slope changes from to 
in Fig. [3]). The lower panel of the same figure shows the 
RDFs obtained at different pressures along the melting 
line. It demonstrates that at low pressure only the steep 
part of the repulsive wall is sampled, whereas at higher 
pressure the softer short-range part of the potential starts 
influencing the atomic interactions. Since softening the 
potential lowers the melting temperatur o^^'^^ we argue 
that the softening of the repulsive interactions in sodium 
at high pressure is the origin of the observed anomaly. 

Our results fit well into the existing theories of met- 
als that predict increase of screening with density and 
possible softening of pair interactions^SiSi. They are also 
in agreement with the fact that certain softening of pair 
interactions can result in anomalous melting^^. 

In order to understand the electronic origin of the soft- 
ening effect we examined the behaviour of the electron lo- 
calization function (ELF), which is widely used to analyse 
the nature of chemical bonds (sec Methods Fig. 2^ 
shows that the behaviour of the derivative of the ELF 
with respect to the interatomic distance, which can be 
taken as a measure of the softness of the bond, closely 
follows the behaviour of the melting curve dramatically 
decreasing at high pressure. It is instructive to compare 
the ELFs at low and high pressures. Fig. shows that 
electrons localized in the Na-Na bond move to the outer 
regions as pressure increases. The tendency of electrons 
to move from the bonding regions to the interstitial re- 
gions at high pressure has already been note d^^i^^ . 

In conclusion, the insight into the electronic and struc- 
tural properties of sodium obtained in this work offer a 
new consistent explanation of the anomalous melting be- 
haviour of Na. We demonstrate that the observed dra- 
matic drop in the melting temperature can be attributed 
to the density dependence of the conduction electrons 
screening, which induces a softening of the interatomic 
potential at short range. These findings have immediate 
implications for explaining behaviour of other metals and 
alloys. 
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FIG. 3. Upper panel: interatomic forces in the region of the 
repulsive wall of the effective pair potential (black solid line). 
Lower panel: RDF of liquid sodium along the melting curve 
at ~10 GPa (blue), ~30 GPa (green) and ~90 GPa (red). 
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FIG. 4. a. Pressure dependence of the derivative of the 
ELF at the center of the Na-Na bond wrt interatomic dis- 
tance taken at the distance of the nearest neighbour in the 
bcc lattice. For comparison, the ma^dmum of Tm is around 
Pe ~ 0.0615A~^. The insets show 2D plots of the same deriva- 
tive in the plane perpendicular to the Na-Na bond for pe 
0.05 A"^ (~20 GPa), 0.07 A"^ (~45 GPa) and 0.09 A"^ 
(~90 GPa). b. Difference between the high-pressure and 
low-pressure values of ELF in the direction perpendicular to 
the Na-Na bond, AELF=ELF(90 GPa)-ELF(25 GPa). 
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Methods. A^A'' potential. The NN potcntiaUi was 
created to reproduce the Perdew-Burke-Ernzerhof (PBE) 
density functional energies of hquid sodium and bcc, fee, 
cI16 crystal phases in the P-T region up to 140 GPa and 
1200 K, as described in our previous paper—. An ultra- 
soft pseudopotential with the 2s and 2p semicore elec- 
trons included explicitly as the valence states was used 
to describe the energetics of the high-pressure structures. 
A large plane- wave cutoff of 100 Ry and a dense mesh of 
fc-points were used for all structures to ensure conver- 
gence of the total energy up to 1.5 meV/atom. The NN 
fitting procedure only introduces a small error (RMSE of 
the independent test sets is 0.91 meV/atom) in addition 
to the numerical (convergence) error of the DPT calcu- 
lations. Thus, the NN-potential closely reproduces the 
ab initio potential energy surface and describes all prop- 
erties and processes in high-pressure high-temperature 
sodium with an accuracy comparable with that of the 
PBE functional. The accuracy of the NN potential for 
numerous structural and dynamical properties of all solid 
and liquid phases in the relevant P-T range has been 
demonstrated in our previous paper—. Supplementary 
Fig. 2 shows that, in addition to these properties, the 
NN potential reproduces a "heat- until- it-melts " melting 
curve obtained from ab initio simulations^. 

Phase diagram. The coexistence lines were determined 
by locating points of equal chemical potential for each 
pair of phases in the P-T plane. This was done in three 
steps (see Ref. [27l ). First, we calculated the Helmholtz 
free energy of each pair of phases by thermodynamic inte- 
gration using the Einstein crystal and the Lennard- Jones 
potential as the reference systems^^. In the next step, 
the chemical potentials were evaluated by integrating the 
free energy as a function of the density. Finally, the co- 
existence lines were traced by integrating the Clausius- 
Clapeyron equation using the predictor-corrector scheme 
of Kofke^. 

In the first step, NN-driven MD simulations were per- 
formed for systems containing several hundred atoms 
(1024 atoms - bcc, cI16, liquid; 868 atoms - fee). 
The temperature was controlled using a coloured-noise 
Langevin thermostat that was tuned to provide the op- 
timum sampling efficiency over all relevant vibrational 
modes^S. The time step was set to 2 fs. The integral in 
the first step was evaluated numerically by the Gauss- 
Legendre quadrature with 10 points. At each value of 
the coupling parameter A, the average value of the in- 
tegrand and its statistical error were obtained from a 
100 ps trajectory. In the second step, the state points 
along the isotherms were obtained from NPT simulations 
(2000 atoms - bcc and cI16, 2048 atoms - fee and liquid) 
governed by the Nose-Hoover equations of motion with 
Langevin noise on the particle and cell velocities^. Av- 
eraging over a 200 ps trajectory was performed for each 
state point. The predictor-corrector algorithm was iter- 
ated until the melting temperature had converged to less 



than 2 K, which required 2-4 iterations of 100 ps each. 

The total simulation time required to model the phase 
diagram amounts to ~50 ns clearly demonstrating the 
computational advantage of the NN approach in com- 
parison with the direct ab initio simulation. 

Jellium model. Two sodium atoms were placed in a 
periodically replicated 17.2 A cubic cell containing free 
electrons with a compensating positive uniform back- 
ground charge. To reproduce the compression effect, the 
number of free electrons was varied so that their density 
is in agreement with that of free 3s electrons in com- 
pressed sodium. To determine the pair potential, the 
ground-state energy of the system was computed for a 
series of fixed interatomic distances with the PBE den- 
sity functional and an ultrasoft semicore pseudopoten- 
tial. A 30 Ry plane-wave cutoff, 1000 K smearing tem- 
perature and a 4 X 4 X 4 Monkhorst-Pack fc-point mesh 
were used to converge the potential curves to 1 mcV. The 
Quantum-Espresso package was used to perform all elec- 
tronic structure calculations. The calculated DFT ener- 
gies were fitted with the analytical form in Eq. [T] The 
fitting parameters are shown in Supplementary Table 1. 

RDFs obtained with the effective pair potential (Sup- 
plementary Fig. 4) perfectly reproduce those from the 
NN-driven simulations. Above ~100 GPa (i.e. in the 
region of the cI16 phase) the pair potential fails to repro- 
duce the more accurate NN simulations because of the 
increasing contribution of many-body interactions. 

Electron localization function. The ELF is calculated 
following the definition given in Ref. [23l 



ELF(r) 



1 



l + {D{r)/Dh{r)Y' 



(2) 



where 



is a measure of the electron localization and 



Du{r)^^^{in'r'p{rf'' 



(4) 



is the reference uniform electron gas value. 

High ELF values show that at the examined position 
the electrons are more localized than in a uniform elec- 
tron gas of the same density, whereas ELF(r) =; 0.5 in- 
dicates that the localization effect is the same as in the 
uniform electron gas. 
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